library(tidyverse)
set.seed(42)
n <- 1000
x <- seq(0, 3, length = n)
Ey <- (sin(3 * x) + 2 * sin(2 * x^2) + 3 * sin(x^3)) / exp(x)
y <- rnorm(n, mean = Ey, sd = 0.2)
mar <- c(2, 2, 2, 2)
xlim<-c(0,3)
ylim <- c(-1.75, 2.5)
cex <- 0.5
par(mar = mar)
plot(x, y, pch = 19, col = "grey", xlim = xlim, ylim = ylim, cex = cex)
lines(x, Ey, lwd = 2)
for(degree in c(5, 10)) {
X = poly(x, degree, raw = TRUE)
plot(x, y, type = "n", xlim = xlim, ylim = ylim,
main = paste0("Polynomial basis functions: degree ", degree))
abline(h = 1, lwd = 2)
for(i in 1:ncol(X)) lines(x, X[,i], col=i+1, lwd=2)
}
for(degree in c(5, 10, 50, 200)){
X = poly(x, degree, raw = TRUE)
model = lm(y ~ X)
plot(x, y, type = "n", ylim = ylim,
main = paste("Polynomial fitted model: degree", degree))
points(x, y, pch = 19, col = "grey", cex = cex)
lines(x, fitted(model), col = "blue", lwd = 2)
legend("topright", paste("AIC =", signif(AIC(model), 3)))
}
for (degree in c(1, 3)) {
n_knots <- 12
df <- n_knots + 1
knots <- seq(0, 3, length = df)[1:n_knots]
X <- 1:n_knots %>%
map_dfc(function(k){
if_else(x < knots[k], 0, (x - knots[k])^degree)
}) %>%
as.matrix
plot(x, y, type = "n", ylim = ylim,
main = paste0(main = "TPS basis functions: degree ", degree, ", knots ", n_knots))
for (k in 1:n_knots) {
lines(x[X[, k] != 0], X[X[, k] != 0, k], col = k + 1)
}
abline(h = 1)
abline(v = knots, lty = 2, col = "grey")
}
for (n_knots in c(12, 24)) {
for (degree in c(1, 3)) {
df <- n_knots + 1
knots <- seq(0, 3, length = df)[1:n_knots]
X <- 1:n_knots %>%
map_dfc(function(k){
if_else(x < knots[k], 0, (x - knots[k])^degree)
}) %>%
as.matrix
model <- lm(y ~ X)
plot(x, y, type = "n", ylim = ylim,
main = paste0(main = "TPS fitted model: degree ", degree, ", knots ", n_knots))
points(x, y, pch = 19, col = "grey", cex = cex)
lines(x, fitted(model), col = "blue", lwd = 2)
legend("topright", paste("AIC =", signif(AIC(model), 3)))
}
}
make_knots = function(degree, n_knots_internal) {
knots <- seq(0, 3, length = n_knots_internal + 2)
knots[-c(1, length(knots))]
}
make_X = function(degree, knots) {
df <- length(knots) + degree + 1
X = splines::bs(x = x, knots = knots, degree = degree, intercept = TRUE)
as.matrix(X)
}
plot_X = function(X, knots = NULL, main = "") {
plot(x, y, type = "n", ylim = ylim,
main = paste0(main = "B-spline basis functions: degree ", degree, ", internal knots ", length(knots)))
for (k in 1:ncol(X)){
lines(x[X[, k] != 0], X[X[, k] != 0, k], col = k + 1)
}
if (!is.null(knots)) {
abline(v = c(0, knots, 3), lty = 2, col = "grey")
}
}
for(n_knots_internal in 11){
for(degree in c(1, 3)){
knots <- make_knots(degree, n_knots_internal)
X <- make_X(degree, knots)
plot_X(X, knots, main = "")
}
}